knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(raster)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
source(here('common_fxns.R'))
library(animation)
reload <- FALSECreate a species richness map for threatened marine species.
Combine taxa-level impact maps to create maps of number of threatened species impacted per cell, per year. This is also done by range priority (see the get_spp_priority() function for details)
Plot maps as animated rasters, both by absolute numbers of species impacted in a location and by proportion of species (impact counts / species richness).
This is all done at the level of cumulative impact category, not individual stressors.
Set up the dataframe of species to include. By taxon, load all species range maps, smash down to taxon-level species richness; then combine these for a total species richness.
spp_incl <- get_incl_spp() %>%
mutate(taxon = tolower(rank))
spp_sens <- spp_incl %>%
filter(!is.na(stressor))
# spp_incl$iucn_sid %>% n_distinct()
# spp_incl %>% filter(!is.na(stressor)) %>% .$iucn_sid %>% n_distinct()
### 1036 if we limit to just sensitive spp, 1271 if we include all
### spp that are mapped/threatened/comp assessedcell_id_rast <- raster(here('_spatial/cell_id_mol.tif'))
ocean_a_rast <- raster(here('_spatial/ocean_area_mol.tif'))
spp_rich_mapfile <- here('_output/rasters/n_spp_map.tif')
if(!file.exists(spp_rich_mapfile) | reload) {
spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
taxa_gps <- spp_incl$taxon %>% unique() %>% sort()
taxa_map_list <- vector('list', length = length(taxa_gps))
for(i in seq_along(taxa_gps)) {
# i <- 2
taxa_gp <- taxa_gps[i]
spp_ids_taxa <- spp_incl %>%
filter(taxon == taxa_gp) %>%
.$iucn_sid %>%
unique()
map_list <- parallel::mclapply(spp_ids_taxa, mc.cores = 12,
FUN = function(id) {
# id <- spp_ids_taxa[1]
spp_f <- file.path(spp_range_dir, sprintf('iucn_sid_%s.csv', id))
spp_df <- read_csv(spp_f, col_types = 'ii') %>%
mutate(iucn_sid = id) %>%
mutate(present = (presence != 5) & !is.na(presence))
})
map_df <- bind_rows(map_list)
map_nspp_df <- map_df %>%
group_by(cell_id) %>%
summarize(n_spp = sum(present))
taxa_map_list[[i]] <- map_to_rast(map_nspp_df, cell_val = 'n_spp')
}
taxa_map_stack <- stack(taxa_map_list)
spp_rich_rast <- raster::calc(taxa_map_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(spp_rich_rast, spp_rich_mapfile, overwrite = TRUE)
}
spp_rich_rast <- raster(spp_rich_mapfile)
plot(log(spp_rich_rast),
axes = F, col = hcl.colors(n = 20),
main = 'log(Species richness)')For each impact category (including all):
calc(fun = sum, na.rm = TRUE)spp_range_dir <- file.path(dir_bd_anx, 'spp_rasts_mol_2020')
impact_map_dir <- here('_output/rasters/impact_maps')
taxon_impact_map_dir <- file.path(dir_bd_anx, 'taxon_impact_rasts')
all_maps <- list.files(taxon_impact_map_dir,
pattern = '.tif',
full.names = TRUE)
taxa_map_df <- data.frame(map = all_maps) %>%
mutate(year = str_extract(basename(map), '[0-9]{4}') %>% as.integer(),
impact = str_replace_all(basename(map), '.+[0-9]{4}_|.tif', ''))make_ramp <- function(n, palette) {
if(n > 101) {
n_even <- round(n + 50, -2)
breaks <- seq(.5, n_even, length.out = 100)
labels <- seq(0, n_even, by = 100)
} else {
n_even <- n
breaks <- seq(min(0.5, n / 100), n, length.out = 100)
labels <- seq(0, n, by = (n / 10))
}
colors <- hcl.colors(length(breaks), palette = palette)
### For breaks and colors, add a different color for
### a value of exactly zero
return(list('colors' = c('grey40', colors),
'breaks' = c(0, breaks),
'labels' = labels))
}
make_gifs <- function(map_stack, filename, layer_names = NULL) {
if(is.null(layer_names)) layer_names = names(map_stack)
ramp <- make_ramp(max(maxValue(map_stack)), palette = 'viridis')
capture.output({
saveGIF({
for(i in 1:nlayers(map_stack)){
plot(map_stack[[i]],
col = ramp$colors,
breaks = ramp$breaks,
axes = FALSE,
axis.args = list(at = ramp$labels),
main = layer_names[i])
}},
interval = 0.5, movie.name = filename,
ani.width = 700, ani.height = 420)
})
return(invisible(NULL))
}
nspp_rast <- spp_rich_rast
values(nspp_rast)[values(nspp_rast) == 0] <- NAstr_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
yrs <- 2003:2013
outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
impact_map_files <- taxa_map_df %>%
filter(impact == str_cat)
for(y in yrs) {
# y <- 2003
outfile <- sprintf(outfile_stem, str_cat, y)
# message('Processing ', basename(outfile))
if(!file.exists(outfile) | reload) {
imp_yr_files <- impact_map_files %>%
filter(year == y) %>%
.$map
impact_yr_stack <- stack(imp_yr_files)
impact_yr_rast <- raster::calc(impact_yr_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(impact_yr_rast, outfile, overwrite = TRUE)
} #else {
# message('File exists: ', outfile, '... skipping!')
# }
}
}outfile <- file.path(impact_map_dir, 'impact_all_2013_2plus.tif')
if(!file.exists(outfile) | reload) {
impact_map_files <- taxa_map_df %>%
filter(impact == 'all_2plus')
impact_2plus_stack <- stack(impact_map_files$map)
impact_2plus_rast <- raster::calc(impact_2plus_stack,
fun = sum, na.rm = TRUE) %>%
mask(ocean_a_rast)
writeRaster(impact_2plus_rast, outfile, overwrite = TRUE)
}for(str_cat in str_cats) {
# str_cat <- str_cats[1]
gif_file <- here(sprintf('figs/impact_pct_all_spp_%s.gif', str_cat))
if(!file.exists(gif_file) | reload) {
rast_files <- sprintf(outfile_stem, str_cat, yrs)
map_stack <- stack(rast_files)
x <- map_stack / nspp_rast
if(any(values(x) > 1, na.rm = TRUE)) {
stop('pct raster values greater than 1, wtf')
}
values(x)[values(x) > 1] <- 1
x <- x * 100
make_gifs(x,
filename = gif_file,
layer_names = sprintf('Impacted spp (percent): %s (%s)', str_cat, yrs))
}
# knitr::include_graphics(gif_file)
cat(sprintf('', gif_file))
}str_cats <- c('fishing', 'climate', 'land-based', 'ocean', 'all')
outfile_stem <- file.path(impact_map_dir, 'impact_%s_%s.tif')
for(str_cat in str_cats) {
# str_cat <- str_cats[1]
r <- raster(sprintf(outfile_stem, str_cat, 2013))
# r_pri <- raster(sprintf(outfile2_stem, str_cat, 2013))
r_norm <- r / nspp_rast
# r_pri_norm <- r_pri / priority_sum_rast
# plot(r, col = hcl.colors(n = 20),
# axes = FALSE,
# main = paste('count spp impacted:', str_cat, 'stressors'))
plot(r_norm, col = hcl.colors(n = 20),
axes = FALSE,
main = paste('pct spp impacted:', str_cat, 'stressors'))
# plot(r_pri_norm, col = hcl.colors(n = 20),
# axes = FALSE,
# main = paste('priority pct spp impacted:', str_cat, 'stressors'))
}r <- raster(sprintf(outfile_stem, 'all', '2013_2plus'))
r_norm <- r / nspp_rast
plot(r_norm, col = hcl.colors(n = 20),
axes = FALSE,
main = paste('pct spp impacted: all stressors (2+)'))for(str_cat in str_cats) {
# str_cat <- str_cats[1]
r_stem <- here('_output/rasters/impact_maps/impact_%s_2013.tif')
imp_r <- raster(sprintf(r_stem, str_cat))
df <- data.frame(impact = values(imp_r),
nspp = values(nspp_rast)) %>%
filter(!is.na(nspp))
df1 <- sample_frac(df, .02)
thingplot <- ggplot(df1, aes(x = nspp, y = impact)) +
geom_point(alpha = .2) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
geom_smooth(method = 'lm', color = 'green', alpha = .5) +
labs(title = sprintf('%s impacts vs spp richness (2 pct sample)', str_cat))
print(thingplot)
}